Long‐term cattle grazing shifts the ecological state of forest soils

Abstract Cattle grazing profoundly affects abiotic and biotic characteristics of ecosystems. While most research has been performed on grasslands, the effect of large managed ungulates on forest ecosystems has largely been neglected. Compared to a baseline seminatural state, we investigated how long‐term cattle grazing of birch forest patches affected the abiotic state and the ecological community (microbes and invertebrates) of the soil subsystem. Grazing strongly modified the soil abiotic environment by increasing phosphorus content, pH, and bulk density, while reducing the C:N ratio. The reduced C:N ratio was strongly associated with a lower microbial biomass, mainly caused by a reduction of fungal biomass. This was linked to a decrease in fungivorous nematode abundance and the nematode channel index, indicating a relative uplift in the importance of the bacterial energy‐channel in the nematode assemblages. Cattle grazing highly modified invertebrate community composition producing distinct assemblages from the seminatural situation. Richness and abundance of microarthropods was consistently reduced by grazing (excepting collembolan richness) and grazing‐associated changes in soil pH, Olsen P, and reduced soil pore volume (bulk density) limiting niche space and refuge from physical disturbance. Anecic earthworm species predominated in grazed patches, but were absent from ungrazed forest, and may benefit from manure inputs, while their deep vertical burrowing behavior protects them from physical disturbance. Perturbation of birch forest habitat by long‐term ungulate grazing profoundly modified soil biodiversity, either directly through increased physical disturbance and manure input or indirectly by modifying soil abiotic conditions. Comparative analyses revealed the ecosystem engineering potential of large ungulate grazers in forest systems through major shifts in the composition and structure of microbial and invertebrate assemblages, including the potential for reduced energy flow through the fungal decomposition pathway. The precise consequences for species trophic interactions and biodiversity–ecosystem function relationships remain to be established, however.


| INTRODUC TI ON
Large ungulate herbivores can profoundly modify ecosystems and the structure of plant-animal assemblages (Stark et al., 2000;Wang et al., 2006;Wardle et al., 2001). Interest in cattle grazing of forest has increased recently due to benefits for nature conservation (Van Uytvanck et al., 2014) and trends toward agroforestry (Röhrig et al., 2020) and forest fire prevention (Ruiz-Mirazo, 2011). Wood pasture and forest grazing by cattle are ancient forms of agriculture that shaped the European landscape from at least the Middle Ages up to the 19th century (Kirby et al., 1995), but which were almost completely abandoned in Northwest Europe due to changes in agricultural practices (Bomanowska & Kiedrzyński, 2011).
Herbivorous animals may shift soil nutrient dynamics from slow-cycling, recalcitrant fungal decomposition pathways to faster nutrient turnover via excreta and bacterial decomposition, hence shortcutting the decomposition pathway (Bardgett & Wardle, 2003).
These changes in nutrient cycling may have cascading effects on the soil food web (Wang et al., 2020), affecting microbes and soil fauna in complex, context-dependent ways. Moreover, an increased level of disturbance (e.g., due to trampling), and resulting soil compaction, may reduce water holding capacity (Houlbrooke & Laurenson, 2013), affect root growth (Unger & Kaspar, 1994), alter microbial activity and biomass (Tan et al., 2005), and negatively affect soil mesofauna (Cole et al., 2008), which may affect nutrient cycling. While direct and indirect effects of ungulate grazing change diversity and abundance of most functional groups in the forest soil, it may also lead to a shift in community composition (Behan-Pelletier, 1999).
Increased nutrient cycling (Siepel, 1996) and physical disturbances (Cole et al., 2008) may cause shifts from long-living K-species to fastreproducing r-species.
Compared to the large body of knowledge about the effects of cattle grazing on pasture soils (e.g., Clapperton et al., 2002;Wang et al., 2020;Yang et al., 2013), far less is known about the impact on forest ecosystems. Most studies focus on the effect on vegetation (Tasker & Bradstock, 2006;Van Uytvanck & Hoffmann, 2009) and forest regeneration (Fortuny et al., 2020) or assess the effects of wild browsers such as deer (Popma & Nadelhoffer, 2020). Moreover, while there is research from boreal (e.g., reindeer; Santalahti et al., 2018;Stark et al., 2000) and tropical forests (e.g., cattle; Stern et al., 2002), the effects of managed ungulates on forest soil diversity and functioning in temperate regions remain largely unexplored.
In this study, we performed a comprehensive investigation of how long-term grazing of birch forest patches shifted the soil physicochemical state and community structure (composition, diversity and abundance of taxonomic or functional groups) of soil organisms (microbes, nematodes, collembola, mites and earthworms). A priori, we hypothesized causal relationships between plant-soil ecosystem components and grazing ( Figure 1a) with specific predictions concerning the different soil organisms. We predicted that the presence of cattle would increase soil nutrient content (nitrogen, phosphorus), leading to a shift from fungal to bacterial dominated soil food webs with concomitant increases and decreases in abundance of bacterivorous and fungivorous nematodes, respectively. Additionally, we predicted that grazing would have a negative effect on soil microarthropods, due to decreased pore space and increased disturbance from trampling, as observed in pastures (Schon et al., 2012). For earthworms, we predicted that faster passaging of organic matter (excreta) would benefit overall earthworm biomass, but that this would depend on differential responses of functional groups defined by their ecological niche (epigeic, surface active or anecic, deep burrowing) and vulnerability to physical disturbance of the soil surface. Finally, we expected that these processes would culminate in distinct assemblages of plants and soil organisms in the disturbed and undisturbed habitat.  (Table S1.1; Figure S1.1). Ungrazed K E Y W O R D S collembola, earthworms, forest grazing, oribatid and mesostigmatid mites, soil chemistry, soil microbes

T A X O N O M Y C L A S S I F I C A T I O N
Agroecology; Community ecology; Functional ecology; Soil ecology semi-natural sites that served as the baseline comparison were ungrazed by cattle or other livestock for at least 70-100 years. All sites were, however, exposed to freely ranging wild herbivorespredominantly roe deer, Capreolus capreolus (Linneaus, 1758)-but their effect on the ecosystem was limited compared to the longterm and concentrated impact of periodically confined cattle herds.
Cattle grazing occurred year-round, but with periodic and haphazard rotation of the livestock out of the woodland according to farmer judgement on the vegetation state to allow recovery. During periods of winter grazing, supplementary feed (hay) was provided at a single feeding station per site. Although we avoided sampling at these locations where cattle feeding was concentrated, this supplementary feeding represented an additional carbon input to the soil ecosystem via cattle excrement. No other systematic differences in management (e.g., forestry practices) occurred at the grazed and ungrazed sites.
F I G U R E 1 (a) Hypothesized causal relations between plant-soil ecosystem components and grazing (note: vegetation-nematode relationship was only tested for phytophagous nematodes). (b) Quantitated relationships between biotic and abiotic components of the plant-soil ecosystem from best subsets of GLMs (AICc). Arrows point in the direction of hypothesized causal relationships (− black/+ gray) with their size proportional to standardized coefficients from the optimal linear model. Dashed lines indicate nonsignificant relationships (p > .05) retained after AICc optimization. For visual purposes, abundance (biomass for earthworms) and not species richness is shown for invertebrates. Vegetation characteristics were not important drivers of other biotic properties and were omitted 2.2 | Data collection

| Sampling design
In each site, we established a 25 m × 25 m sampling grid with 16 potential nodes at the center of the woodland (≥50 m from the forest edge). Sampling took place on randomly selected nodes of this grid (4-6 sampling points depending on the variable, with a minimum 5-m separation) to capture within-site spatial heterogeneity. The timing of sampling within the season was designed to coincide with the peak of activity or diversity of each taxon, while the year of sampling (2007 or 2008) was dictated by logistical constraints (see below).
Because sampling points were pseudoreplicates, we pooled data for statistical analysis at the level of the forest patch site (see below).

| Forest understory plant community
The species composition and percentage (%) cover of the herbaceous understory plant community (vascular, bryophytes, and macrolichens) of each site was measured in July 2008 in a series of five quadrats (1 m 2 ) randomly situated on the grid nodes. The identity and visually estimated percentage cover of each species was determined, enabling the plant species richness and percentage cover of broad plant groups (dicotyledonous herbs, graminoids, pteridophytes, and bryophytes) to be derived.

| Soil properties
In 2008, a bulked sample of six soil cores (50 mm diameter, 50 mm depth), minus the litter layer, were taken per site and analyzed for pH (measured in 0.01 M CaCl 2 solution), phosphate levels (colorimetric PO4-P Olsen mg/kg), and the total percentage (%) content of carbon and nitrogen using an Elementar Vario EL elemental analyzer (n = 20 for all). Soil bulk density (g/cm 3 ) was estimated by driving a steel circular ring of standard volume (77 cm 3 ) into the soil (250 mm depth) at three random grid points per site. Samples were oven dried and dry mass of soil and parent material fraction (g −1 ) obtained and converted to a volumetric mass of dry soil per sampling point. We then obtained a mean value per site of all these soil property values for subsequent analysis.

| Microbial community structure
Standard phospholipid fatty acid (PLFA) analysis was used to quantify the dry weight-based mass of markers of fungal and bacterial biomass and the fungal-to-bacterial ratio in the soil samples from the sites, see Frostegård et al. (2011) for a discussion of this method.

| Nematode abundance, richness, and trophic group
Four soil cores (25 mm diameter, depth 50 mm) were taken (May 2008) per site (80 cores in total) and weighed to determine the wet mass (g) per sample. Nematodes were extracted (48 h) using a modified Baermann funnel extraction method (van Bezooijen, 2006). The soil sample was placed into a sieve (110 mm diameter, 1 mm mesh) lined with a filter tissue (KimTech Delicate task wipe, Kimberly-Clark™) and situated within a water-filled funnel (150 ml) to encourage nematode migration from the soil into the water column where they gravitate into a vial (4 ml) connected to the funnel by rubber tubing. van Bezooijen, 2006) followed by desiccation (40°C). We pipetted and slide mounted (wax sealed) a drop of this nematode suspension from each sample. We determined the relative abundance of nematodes to the highest possible taxonomic resolution (Family/Genus) in a community subsample of 100 individuals (×1000 magnification, Leica DM12.5 microscope) and assigned them to a trophic group (e.g., predator, omnivore, bacterivore, or herbivore; Bongers, 1994).
For each subsample, the number of nematodes per trophic group was divided by the dry mass (g) of the individual soil sample and multiplied by 100 to get the number of individuals per 100 g. The values for the four subsamples were then averaged to obtain a mean for each site.
Observed species and their abundance are given in Table S3.1.
To describe nematode community structure, we calculated the Enrichment Index (EI), Structure Index (SI), and Channel Index (CI) using the assignment to nematode trophic groups and colonizer-persister classes (Bongers, 1990). The EI is an indicator of soil nutrient enrichment, the SI indicates food web stability (high SI = stable food webs dominated by disturbance-sensitive species with high longevity; Ferris et al., 2001). The CI indicates the relative importance of the fungal energy channel compared to the bacterial energy channel in the soil food web (Cesarz et al., 2015). Further description of these indices is in S4.

| Earthworms
We sampled earthworms (Lumbricidae) in May 2008 using an iron frame driven into the soil to delineate a monolith of soil (250 × 250 × 150 mm depth). These monoliths were taken from three random points on the sampling grid within each site (n = 60).
After removing the litter layer, we excavated and sorted the soil by hand to a depth of 150 mm, earthworms collected were immediately preserved by immersion in a vessel containing 30 ml 80% Industrial Methylated Spirits solution. The species identity, ecological niche (epigeic/endogeic/anecic), total abundance, and total pooled fresh biomass (g) of earthworms per monolith sample were determined (Sims & Gerard, 1985) and are displayed in Table S3.2.

| Multiple impacts of grazing on the forest plant-soil ecosystem
The effect of the cattle grazing on the birch ecosystem was analyzed using general linear and best-subset modeling (AICc for small sample sizes; Burnham & Anderson, 2002) implemented with the MuMIn package (Barton, 2020). Given the likelihood of multiple impacts of the grazing cattle, the many intercorrelated variables (Table S1. Table S1.2) to assess the direct correlation. Then, we ran GLMs for each category of ecosystem variable (e.g., soil properties; vegetation; microbes; etc.). For this, we ran a two-tier approach, first running full models for each variable category (e.g., soil model, microbial model,…) and selecting the best subset of variables by AICc optimization for each of these categories. Secondly, we compared the model of grazing presence (0/1), all the best models per category, and an intercept-only (null) model to select the final best model based on AICc scores (Tables 1   and S5.1). Models with a ΔAICc < 2 compared to the optimal model were also considered plausible. Highly correlated explanatory variables (e.g., Figure S6.1% C & N) were not fitted to the same models.
We tested for homoskedasticity (Levene) and normality of residuals (Shapiro-Wilks) and, where necessary, the response variables were log-transformed to meet model assumptions.

| Grazing impacts on composition of plant and soil invertebrate communities
We tested how grazing affected the community composition of plants, nematodes, collembolans, oribatid mites, mesostigmatid mites, and earthworms. We used nonmetric multidimensional scaling (NMDS) to visualize differences in community composition based on a Bray-Curtis dissimilarity matrix (500 iterations; R vegan package: metaMDS function Oksanen et al., 2019). Stress scores were sufficiently low to interpret the plots in two dimensions. Data were log-transformed to reduce the influence of the most dominant taxa.
We used a Permutational Analysis of Variance (PERMANOVA, 9999 permutations; R vegan package, adonis2 function) to test if there was a significant difference in community composition between the grazed and ungrazed sites with homogeneity of variance tested using the betadisper function.

| Forest soil properties
Cattle grazing led to a nearly four-fold increase in available soil phosphorus levels (Table 1; Figures 1b and 3a). While no significant differences in the percentage of soil carbon and nitrogen were detected between grazed and ungrazed forests, the C:N ratio was significantly lower in grazed sites (

| Forest understory vegetation
Grazing significantly increased herb layer cover and plant species richness (Table 1; Figure 2) in the forest patches. Plant species richness was related positively (Table 1) (Table 1).

| Soil microbial PLFA
The concentration of microbial PLFA was lower in grazed forests, with bacterial PLFA decreasing by almost a third and fungal PLFA being almost 60% lower (.2). This was correlated with the lower C:N ratio in grazed sites (Table 1; Figure 1b), which, compared to the grazing model, better explained the lower bacterial and fungal PLFA (ΔAICc = −3.9 and −6.6 respectively, Table S5.1). The stronger decline in fungal PLFA shifted the ratio of fungal to bacterial PLFA strongly toward bacteria, with bacteria representing 79.1% of all F I G U R E 2 Standardized coefficients from a GLM with grazing treatment as the single explanatory variable of the response of plant-soil ecosystem components to the presence or absence of cattle grazing in forest patches ( NS p > .10; o p < .10; *p < .05; **p < .01; ***p < .001); Veg., Vegetation. Positive values indicate a positive correlation with presence of cattle grazing PLFA in ungrazed plots and 83.6% in grazed plots (Figure 2). The difference in fungal-to-bacterial PLFA ratio was best explained by the "grazing" model (Table 1).

| Nematode trophic groups and food webs
While nematode taxonomic richness marginally increased under Note: Variables that were best explained by an intercept only model (%C, %N, phytophagous nematodes, bacterivore nematodes, omnivorous nematodes, SI, collembolan richness, and endogeic earthworm biomass) are not displayed here; see Table S5.1.

TA B L E 1 (Continued)
F I G U R E 3 Box plots showing the effects of grazing treatment on forest soil properties. Significance levels inferred from linear models using only grazing treatment as a response variable are indicated above each figure (NS: not significant, ¤ p < .10, *p < .05, **p < .01, ***p < .001)  Table 1).

F I G U R E 4 NMDS
Fungivorous nematode abundance was significantly lower in grazed sites with a 75% decrease in numbers compared to the ungrazed sites (Table 1; Figures 1b and 2). Other nematode trophic groups did not show a significant response to grazing (Figure 3) with differences between sites best explained by an intercept-only (null) model, except for predatory nematodes, which showed a positive correlation with PLFA markers of fungal biomass (ΔAICc = −5.0, compared to intercept-only (null) mode, Tables 1 and S5.1).
Grazed forests had a significantly higher EI and a significantly lower Channel Index, while the SI did not differ significantly between grazed or ungrazed forests (Table 1; Figure 2).

| Microarthropod abundance and richness
The presence of grazing cattle had a strong negative direct influence on the abundance and species richness of oribatid and mesostigmatid mites (Figure 2). Collembolan abundance was associated negatively with the presence of grazing, but species richness was unaffected (Table 1; Figure 2).
Oribatid abundance was best explained by a soil-only model, containing Olsen phosphorus, soil pH and soil bulk density (Table 1). This model was considered equivalent to the grazing model

| Earthworm biomass
Total earthworm biomass was significantly higher in grazed forest patches (Table 1; Figures 1b, 3 and S6.2). However, this effect of grazing on earthworm biomass varied among functional groups.
Anecic earthworms were completely absent from ungrazed plots, but the most dominant guild in grazed forests, where they represented 42.5% of earthworm biomass (Table 1; Figure 5). While grazing best explained this pattern, a soil model containing the C:N ratio was also considered (ΔAICc = −1.6). Epigeic and endogeic earthworms were unaffected by grazing ( Figure 5). Epigeic earthworm biomass was positively and (marginally) negatively correlated with PLFA markers of bacterial and fungal biomass, respectively (Table 1)

| Soil invertebrate assemblage composition
NMDS and PERMANOVA analysis revealed clearly distinct assemblages of nematodes, collembolans, oribatids, mesostigmatids, and earthworms associated with the grazed or ungrazed forest patches ( Figure 4b-f). While the oribatid assemblages showed a much lower dispersion in ungrazed plots (p = .0001), separation between the two groups was sufficiently strong to consider the results reliable, given the robustness of this test for balanced designs (Anderson & Walsh, 2013).

| DISCUSS ION
Long-term cattle grazing consistently transformed the plant-soil assemblages in these replicated forest patches by modifying the abiotic soil environment and directly-or indirectly via changes in the soil physicochemical state or biological associations-the abundance, diversity, and composition of soil biota. One major functional consequence of the long-term introduction of cattle to this ecosystem was the reduction of the importance of the fungal energy channel, lower microarthropod abundance and diversity, and a functional shift in hence litter inputs (Bardgett & Wardle, 2003;Grayston et al., 1996;Hamilton et al., 2008).
The soil physical structure was similarly modified in the grazed habitat with slightly elevated soil bulk density implying compaction from trampling by herds, which reduces soil porosity and water retention capacity (Sharrow, 2007).
Carbon is the limiting nutrient for microbial detritivores outside the rhizosphere (Grayston et al., 1996;Scheu et al., 2005). Here, cattle grazing reduced the overall soil microbial biomass, which corresponded with the fall in the C:N ratio. Similar declines in microbial biomass have been reported in European submontane grassland pastures under intensive grazing (Bardgett et al., 2001), although in subtropical pastures, grazing can increase microbial biomass (Wang et al., 2006), pointing to how environmental zone can introduce context dependencies. Additionally, since microbial biomass might respond in a nonlinear way, a gradient in grazing intensity might reveal patterns obscured by our binary approach (Bardgett et al., 2001).
Although N-saturated soils may show declines in soil microbes (Wallenstein et al., 2006), we found no evidence that grazing shifted the soil N level.
Fungal biomass was more negatively affected than the bacterial component, potentially because fungi are more associated with soils with high C:N ratios (Oates et al., 2012;Wan et al., 2015), recalcitrant litter (Högberg et al., 2007), and mycelia are more susceptible to damage from physical disturbances, that is, trampling and soil compaction (Hartmann et al., 2014). The shift in the dominance structure of the soil microbial community, here caused by a reduction in fungal biomass in the grazed sites, is consistent with previous studies that show grazing-induced shifts toward the bacterial energy channel (Lopez-Sangil et al., 2011;Oates et al., 2012;Waring et al., 2013). Shifts from fungal to bacterial decomposition have been associated with changes in carbon sequestration (Malik et al., 2016), but the magnitude of the effect in this study system was insufficient to alter soil carbon content.
This shift to more bacteria-dominated soil food webs in grazed forest patches mediated through changes in nutrient ratios led to corresponding knock-on effects on the soil nematode assemblage.
Fungivorous taxa were the only nematode trophic group to be significantly less abundant in grazed forest patches, but a significant positive correlation between fungal biomass and the presence of predatory nematodes suggests a possible bottom-up conduit in the nematode food web connecting fungivorous nematodes and their predators (Wardle & Yeates, 1993). The significantly lower nematode Channel Index in grazed forest fragments indicates a shift from fungivorous to bacterivorous nematode dominance. This provides further evidence of a decreased energy flow through the fungal channel in the grazed forest patches. This is in line with studies on fertilized grasslands (Parfitt et al., 2010) and links to relatively higher soil N content (Cesarz et al., 2015), but contrasts with other studies that found a positive correlation, or no correlation at all between grazing and the Channel Index in grassland ecosystems (Briar et al., 2012;Mills & Adl, 2011). This is likely due to the fact that the bacterial energy channel usually predominates in intensively grazed grasslands, while forest soils generally have a dominant fungal energy channel (Ruess, 2003). While different trophic groups contrasted in their responses to grazing and linked environmental characteristics, total nematode richness was positively correlated with relatively greater soil pH in grazed habitat (>5.0) and consistent with the optimal pH range (5.0-7.0) for most nematode genera (Mulder et al., 2003;Nisa et al., 2021).
This shift in the dominance of nematode groups was also reflected in a grazing-associated increase in the EI, caused by faster nutrient cycling through input of more labile organic matter and a proliferation of more opportunistic bacterial-feeding nematodes (Ferris et al., 2001). This response of the nematode assemblage structure to enrichment can be idiosyncratic with positive and negative responses seen in grassland systems often dependent on grassland habitat type and grazing intensity (Hu et al., 2015;Wang et al., 2006Wang et al., , 2018. The observed nutrient-driven shifts in nematode assemblages to a bacterial energy channel did not translate, however, into an overall reduction in the stability or complexity of the nematode food web (SI; Ferris et al., 2001). Trophic linkages were rarely detected in our models. This may reflect the spatial scale at which we undertook our analysis, with effects attributable to grazing detectable, but lacking the finer spatio-temporal resolution to detect the signal of soil species interactions and their cascading effects across the soil food web (Thakur et al., 2020). Consequently, subtler trophic interactions went undetected and only the modified nutrient cycling that shifted the relative dominance of nematode feeding groups in the microfood web was observed.
Microarthropod abundance and species richness (excepting Collembola) were strongly reduced by grazing in these forest patches, which was linked to changes in the soil physicochemical state (elevated phosphorous content, bulk density and pH). By inhabiting the interstitial soil pore space, microarthropods may have suffered greater direct mortality from soil compaction by cattle trampling as well as loss of microniche space or nutritional resources with the physical disturbance of the litter layer (Bardgett et al., 1998;Hopkin, 1997;Wardle et al., 2001). Larger bodied arthropods are more sensitive to the adverse effects of grazing (and other disturbances) than small-bodied species (Wardle & Bardgett, 2004). Here, due to their small size nematodes were unaffected by soil compaction (Bouwman & Arts, 2000;Schon et al., 2012), while mites showed a far stronger negative response than collembolans (Andriuzzi & Wall, 2017).
In particular, oribatid mites are known to be associated with acidic soils (Maraun & Scheu, 2000) and are very sensitive to physical disturbances (Schon et al., 2012) with low resilience due to longer generation times and "K-selected" traits (Cole et al., 2008). Collembolans were also more abundant on the more acidic soils of the ungrazed habitat, possibly due to physiological adaptation to that niche or their role as grazers of the fungi typically associated with more acidic soils such as under forest (Hopkin, 1997;Vanbergen et al., 2007), although we found no direct correlation between collembolan abundance and PFLA markers of fungal biomass here.
Grazing also produced distinct microarthropod assemblages.
Mites, particularly oribatids, and to a lesser extent collembolans, showed a much higher dispersion of community composition in grazed sites compared to the relatively homogenous communities in ungrazed forest patches. This suggests that cattle grazing in these forest patches created a higher level of beta-diversity between sites. The mechanism behind this remains to be established.
However, it would be consistent with greater soil microhabitat or niche diversity arising from the actions of cattle on the plant-soil ecosystem.
As predicted, the response of earthworms to grazing was contingent on their exact ecological role and function. Epigeic and endogeic earthworm biomass did not differ significantly between grazed and ungrazed sites, this was unexpected for epigeic species that are vulnerable to the effects of trampling .
Manure can be a food resource for earthworms (Bacher et al., 2018;Curry et al., 2008) and hence the addition of this partly digested organic matter input to these forest soils may have compensated for adverse effects associated with any physical disturbance by the grazing livestock. Strikingly, anecic earthworms were absent from ungrazed sites, whereas they were the most dominant earthworm functional group in grazed forest patches. Because of their vertical burrowing behavior, anecic earthworms can benefit from manure input at the soil surface, while their deep vertical burrows protect them from trampling by cattle (Schon et al., 2012). This may explain the observed increase in anecic earthworm biomass associated with grazing, as reported from grassland systems (Muldowney et al., 2003;. This has implications for ecosystem functioning as anecic earthworms vertically transport nutrients and organic material (Don et al., 2008) and increase nitrogen mineralization (Van Groenigen et al., 2014). Although they may compensate for soil compaction caused by cattle (Capowiez et al., 2012) and form more stable microhabitats for soil arthropods, rich in nutrients and microorganisms (Eisenhauer, 2010), in this system, their uplift in biomass was insufficient to offset the other impacts.
Overall, habitat disturbance of birch forest patches by long-term grazing by ungulates profoundly modified the plant-soil subsystem probably via the combination of trampling, biomass consumption, and excreta. Specific shifts in nematode community trophic structure indicated a reduction in the predominance of energy flow through the fungal pathway and a subtle shift toward a soil food web underpinned by bacterial decomposition of more labile organic matter. Such results demonstrate the ecosystem engineering potential of large ungulate grazers in forest systems, however, the consequences for biodiversity-ecosystem function relationships (e.g., nutrient cycles, soil-atmosphere gas exchanges) and species interactions remains to be understood.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
Raw ecological data and R-code: Dryad https://doi.org/10.5061/ dryad.wm37p vmq2.